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Abstract. A current-carrying resonant nanoscale device, simulated by non-adiabatic 
molecular dynamics, exhibits sharp activation of non-conservative current-induced forces with 
bias. The result, above the critical bias, is generalized rotational atomic motion with a large 
gain in kinetic energy. The activation exploits sharp features in the electronic structure, and 
constitutes, in effect, an ignition key for atomic-scale motors. A controlling factor for the effect 
is the non-equilibrium dynamical response matrix for small-amplitude atomic motion under 
current. This matrix can be found from the steady-state electronic structure by a simpler static 
calculation, providing a way to detect the likely appearance, or otherwise, of non-conservative 
dynamics, in advance of real-time modelling. 

Nanoscale conductors [1,2] carry current densities orders of magnitude larger than in 
macroscopic wires, resulting in substantial forces. The current-induced force on a nucleus 
consists of (I) the average force, and (II) force noise originating from the corpuscular nature 
of electrons [3,4]. II causes inelastic electron-phonon scattering and Joule heating [2,5]. I 
contains the so-called electron- wind force, and velocity-dependent forces [3,4,6-11]. The 
wind force results from momentum transfer in elastic electron-nuclear scattering, and drives 
electromigration [12]. Wind forces can be calculated from first principles, to study their effect 
on nanoscale devices [13-16]. 

The electron-wind force is receiving fresh attention due to a remarkable property: it is 
non-conservative (NC) and can do net work on atoms around closed paths [3,4, 6-12]. This 
phenomenon, which we call the waterwheel effect, opens up interesting questions. It provides 
a mechanism for driving molecular engines [17-22]. But the gain in kinetic energy of the 
atoms from the work done by NC forces is also a potential failure mechanism, possibly 
more potent than Joule heating. There are experimental indications of anomalous bias- 
activated apparent heating in point contacts [23, 24], above that expected from Joule heating 
alone [13,25-27]. The waterwheel effect is a possible activation mechanism also for the 
electromigration phenomena that become a central issue under large currents [28, 29]. 

The applied bias is a key factor for the operation of NC forces [3,4,6,8-11]. First, 
they have to compete against the electronic friction (a velocity-dependent force), and this 
may require a critical current. Second, the waterwheel effect requires pairs of normal modes 
degenerate in frequency [3,6,8]. If there is a frequency mismatch, a critical bias may be 
needed to overcome it. Ramping up the bias to overcome these factors is, notionally, like 
having to press the accelerator harder to climb a hill. 
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The bias is a sensitive control parameter for the current-driven excitation of atomic 
motion in resonant systems [30, 31]. In this letter we show how resonances can be exploited 
to turn the NC force on and off, akin to turning the engine of a car on and off. We will see that 
this "switch" is robust against decreasing resonance width. We will illustrate further how NC 
forces can be gauged to extract work from the current with no net angular momentum transfer 
to the real-space atomic motion. 

Our system is shown in Fig. 1 . It is a 2D metallic nanowire in the x-y plane. The wire 
has a simple square lattice structure with a bondlength of 2.6 A. Electrons are described in 
a nearest-neighbour single-orbital tight-binding model with parameters for gold [32], except 
the band filling which here is set to 0.16. The hopping integral is H = —3.3175 eV and its 
derivative with distance is H' = 5.1038 eV/A. Red denotes metallic atoms, with an onsite 
energy set to zero. Region C is a resonant device created by the blue atoms, which have an 
elevated onsite energy, E^. These insulating atoms form a double constriction. Current is 
supplied by the leads L and R. 

First we examine the steady-state transport properties of the system, in the Landauer 
picture summarized in Fig. 1. Under bias V, the steady- state 1 -electron density matrix can be 
written as 

/+oo 
p(E)dE, (1) 
-oo 

where p(E) = f L (E)D L (E)+f R (E)D R (E) and D L ^(E) are the density of states operators 
for the stationary scattering states {ipum}, an d R denotes ionic coordinates. We consider 
non-interacting electrons throughout. Spin is subsumed in Dnm{E). 




Figure 1. The 2D atomic strip considered in our calculations. A device C is connected 
to two semi-infinite electrodes L and R. Red denotes metallic regions, the blue atoms are 
insulating, and black denotes the atoms that will be treated dynamically in the simulations. In 
the Landauer picture, particle reservoirs inject electrons from the left and from the right with 
electrochemical potentials = + eV/2 and /j,r = fJ, — eV/2 respectively, where \i is the 
equilibrium Fermi level, e(< 0) is the electron charge and V is the applied bias. eV > 
corresponds to electron flow from left to right. Electrons in the steady state are described by 
Lippmann-Schwinger scattering states, {4>l{r) }, populated with the Fermi-Dirac distributions 
Il(r) corresponding to ^i(ii)- 
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Figure 2. Curl of the current-induced force on atom (1) from Fig. 1 vs Fermi level, to lowest 
order in the bias V. Under finite bias, the curl is given by the integral of the curve (with the 
integrand now in units of 2/ A 2 ) from to /i^. The inset shows the transmission function 
for the system. For the parameters in the text, the equilibrium Fermi level is denoted by the 
black line. The energy window for conduction under several biases (used in the simulations) 
is marked by coloured lines. The onsite energy shift on the barrier atoms is E\, = 10.8 eV 
(red) and E h = 20 eV (blue). 



Forces exerted by electrons on ions are described within the Ehrenfest approxima- 
tion [33]. In the steady state, this force is given by 

F(V,R) = Tr{P(R)p(V,R)}, (2) 

where F(R) = —VRH e (R) is the force operator, defined in terms of the electronic 
Hamiltonian H e (R). The curl of the force on an ion with position (R x , R y , R z ), 

8F BF 

(V x F) 2 = ^ - ^ , (3) 
v ; dR x dR y ' 

is given by [6, 7] 

(VxF) 2 = 4tt f ImTr{F x D(E)F y p(E)\ dE . (4) 

J — oo 

Here, D(E) = D L (E) + D R (E) = [G~(E) - G + (E)]/2ni is the total density of states 
operator, where Ct (-E 1 ) are the retarded and advanced Green operators. The curl comes solely 
from the non-equilibrium part of the density matrix. Thus, at zero electronic temperature 

(V x F) 2 = 4vr Im Tr ^F x D(E)F y D L (E)^ dE 

- 4tt r Im Tr {P x D(E)F y D R (E)\ dE . (5) 

To lowest order in the bias, this simplifies to 

(V x F) 2 = 47rImTr{itDGu)F y Ap} , (6) 

where Ap = eV[D L (fi) - D R (fi)]/2. 

Fig. 2 shows the energy-resolved curl of the force on atom (1) from Fig. 1 and the 
transmission function for the system (inset), for two values of E h . By symmetry, the curl 
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for atom (3) is identical. To see this, reverse the bias. Eq. (6) changes sign. Atom (1) under 
the new bias is physically equivalent to atom (3) under the old bias. It follows that the curl on 
(3) and (1) is the same, under a given bias. The curl on (4) is the negative of that on (1), and 
so on. 

At the given Fermi level, the perfect 4-atom wide strip has two open conduction channels, 
corresponding to the two lowest transverse modes in the wire. The double constriction 
formed by the insulating atoms constitutes a tunnelling double barrier for the higher mode. 
The resultant electronic resonance is present both in the curl and in the transmission in 
Fig. 2. Raising E h makes the walls harder and the effective constriction narrower, making 
the resonance narrower. The shift is the energy renormalization that accompanies changes in 
lifetime for quasibound states. (The effect eventually saturates with E h .) 

First we study the system with the wider resonance, for = 10.8 eV. The peak in the 
curl is at an energy of — 2.21 \H |. Then we expect that we can switch the NC force on and off 
by varying the bias so as to just include, or exclude, the peak from the conduction window. 
For the present system the equilibrium Fermi level is /i — — 2. 1 1 18| H | . The peak falls in the 
voltage range —1.0 V < V < —0.3 V. We thus expect a threshold bias somewhere in this 
range, at which the waterwheel effect is activated. 

To test these ideas we simulate this system numerically. As in Ref. [6], the simulations 
are carried out in the Ehrenfest approximation. Ehrenfest dynamics consists of the quantum 
Liouville equation for the electronic density matrix, coupled with Newtonian equations for 
nuclei. It includes the mean forces (I), but excludes the noise (II) [4,34]. The resultant 
suppression of Joule heating enables us to attribute any observed heating to work done by 
the NC forces. Current is generated by the multiple-probes open-boundary method [35], with 
parameters: S is a 4 x 151-atom strip; C contains 28 atoms; the source and sink terms are 
applied to the 248 end atoms in each electrode; T = 0.6 eV, A = 0.001 eV. 

First we allow only atom (1) from Fig. 1 to move, setting its mass to M — 1 atomic mass 
unit (amu). Fig. 3 shows its kinetic energy in time for a bias of —1.0 V. Although the bias 
window encompasses the peak in the curl of the current-induced force, no heating occurs. 
Instead cooling under the electronic friction is observed. Velocity-dependent forces can be 
calculated perturbatively under steady-state conditions [3,4,9-11]. Another way of isolating 
their effect is via the atomic mass [6]. Increasing M suppresses them relative to the velocity- 
independent NC force, enabling the latter to dominate. This is illustrated in Fig. 3 for M = 20 
amu. 

Next we allow atoms (1) and (3) to move, with M — 1 amu. Fig. 4 presents the total 
kinetic energy of the two atoms for biases of —0.7 V, —0.8 V and —1.0 V. For —0.7 V we see 
only cooling. At —0.8 V the kinetic energy starts to grow, showing that the waterwheel effect 
has been activated. Increasing the bias to — 1 .0 V shows an even larger energy gain. If we 
examine the trajectories we observe unidirectional rotational orbits that spiral outward as the 
energy grows. These are similar to those in bent atomic chains [6], although here they tended 
to be elliptical. 

Thus for the two moving atoms we see the expected bias threshold for the waterwheel 
effect, between —0.7 V and —0.8 V. But for the same mass, for one moving atom the effect 
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Figure 3. Kinetic energy of atom (1) from Fig. 1 during current flow over a period of 3000 fs, 
under a bias of —1.0 V, for the wide resonance (E\, = 10.8 eV). In this simulation only atom 
(1) is allowed to move. The energy window for conduction encloses the peak in the curl of 
the current-induced force in Fig. 2 but when the mass of the atom is M = 1 amu (blue) non- 
conservative heating does not occur. However, when M = 20 amu (red) the non-conservative 
force overcomes the electronic friction and the waterwheel effect is seen. 



_ V = -0.7V E b = 10.8 eV 
_ V = -0.8V E b = 10.8 eV 




Time (fs) 

Figure 4. Total kinetic energy of atoms (1) and (3) in Fig. 1 for the wide resonance (Eb = 
10.8 eV) during current flow over 1000 fs, for three applied biases: —0.7 V, —0.8 V (main 
panel) and —1.0 V (inset). When the magnitude of the bias is below 0.8 V we see only 
cooling, whereas for larger biases the heating due to non-conservative forces dominates. The 
inset shows the kinetic energy for the two atoms at —1.0 V for both the wide resonance and 
the narrow resonance (£t> = 20 eV) and we see that the kinetic energy gain for the narrow 
resonance exceeds that for the wide resonance. M = 1 amu. 



fails to manifest itself even above this threshold. Both features can be understood in terms of 
the non-equilibrium dynamical response matrix. 

We write the dynamical response matrix under current as 

^>*> = -|£ + aSb (7) 

where F v is the electronic force of Eq. (2), Vjj is the ion-ion interaction, and index v labels 
ionic degrees of freedom. K vv > may be split into an equilibrium part K^, and a non- 
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equilibrium correction AK w i, K„ U /(V, R) = K^, + AK VV >. Since V x F is non-zero under 

current, Eq. (3) implies that AK uv i is in general not symmetric and hence mode frequencies 

under current will in general be complex. Complex frequencies imply motion that grows or 

decays in time. The growing modes are the waterwheel, or runaway [3], modes. 

Differentiating the steady-state density matrix by scattering theory as in [6, 7], we obtain 

AK vu t = S vu t + A vv i, where the symmetric and antisymmetric components are given by 

AK VV , + AK V , V 
i>u* - 2 

= Yl r Tr{iWA(£)} dE 

i=L,R P 

+ 2 J2 Re Tr {FuR(E)F u , Di(E)\ dE (8) 

AK VV , - AK V , V 
' ~ 2 

= 27r J2 ImTrjF^^^A^)} dE . (9) 

Above, R(E) = [G~(E) + G + (E)}/2 and K vv , = d 2 H e /dR u dR u < = -dF u /dR u <. 

The competition between S and A will determine whether mode frequencies are real or 
complex. We have calculated these frequencies for the systems above. When only atom (1) 




Figure 5. Total kinetic energy (top) and individual trajectories (bottom), with atoms (l)-(6) 
allowed to move, for the wide resonance (Eh = 10.8 eV). The bias is —1.0 V and M = 20 
amu. The trajectories cover an advanced stage of the motion, with large displacements: this 
interval is denoted in blue in the kinetic energy plot. Dots (squares) denote the start (end) of 
each trajectory. 
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Normal Applied Bias 

frequency Equilibrium -0.7V -0.8V -1.0V 

(fs- 1 ) W.R. N.R. W.R. W.R. W.R. N.R. 

^ 0.442 0.432 0.368 0.470 + 0.016i 0.477 + 0.014i 0.476 + 0.013i 

lo 2 0.457 0.450 0.464 0.470 - 0.016i 0.477 - 0.014i 0.476 - 0.013i 

u 3 0.473 0.471 0.474 0.482 0.483 0.481 

oo 4 Q.494 0.492 0.481 0.490 0.500 0.504 

Table 1. Mode frequencies at equilibrium and at different applied biases when both atoms (1) 
and (3) in Fig. 1 are allowed to move. Results are shown for both the wide resonance (labelled 
W.R.; E b = 10.8 eV) and for the narrow resonance (labelled N.R.; E h = 20 eV). M is 1 amu. 



in Fig. 1 is allowed to move the equilibrium frequencies for M = 1 amu are u)\ = 0.456 fs 
and cu 2 = 0.477 fir 1 , while for a bias of -1.0 V w 1)2 = 0.480 ± 0.005i fs -1 . 

Table 1 shows the mode frequencies when both atoms (1) and (3) are allowed to move. 
All frequencies are real at —0.7 V, while two frequencies become complex at —0.8 V. The 
waterwheel effect should activate in between, as observed in Fig. 4. The imaginary parts in 
Table 1 at — 1 .0 V are considerably larger than for one moving atom, explaining the greater 
ease with which waterwheel motion occurs. Remarkably, the frequencies for the narrow and 
wide resonance in Table 1 have comparable imaginary parts. We can see this from Fig. 2: 
the peak in the transmission is bounded, but the curl grows in height with decreasing width, 
keeping the area under the curve roughly constant. The reason is evident from Eq. (6): the 
imaginary part comes from Ap and is related to the current; but present also is the density of 
states, which is not bounded. We conclude that, at least within limits, the "switch" for the NC 
force can be made more abrupt by making the resonance sharper. 

The electronic friction should drop as the edges of the bias window move away from 
the resonance [36]. This is reflected in the results for the two moving atoms in Fig. 4. The 
imaginary parts of the frequencies at —0.8 V and —1.0 V (for both the wide and narrow 
resonance) are comparable. But the rate of energy gain increases as \xl and jiR retreat from 
the peak. 

In Table 1 the large variation in u>i around —0.7 V arises from the behaviour of the 
integrals in Eqs. (8) and (9), as the resonance moves into the conduction window. The 
interplay between the components of the dynamical response matrix will be studied further in 
future work. 

Allowing atoms (1), (3), (4) and (6) to move produces only real frequencies at —1.0 
V, but with all six atoms allowed to move two waterwheel modes appear, one with a larger 
imaginary part than any above. We simulate the system, starting with small random velocities, 
in Fig. 5. To slow down the motion, we set M = 20 amu. The trajectories show that there 
is angular momentum transfer to each atom, but no net angular momentum transfer overall. 
Thus, although the elemental motion under NC forces is a generalized rotation in the abstract 
space of two normal coordinates coupled by current [8], these forces can do work without 
angular momentum transfer to the real-space dynamics of the atomic subsystem as a whole. 
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We have simulated dynamically a resonant device and have shown how the bias can 
switch on and off the non-conservative forces on atoms induced by current. This switch 
appears robust with decreasing resonance width. This robustness will ultimately be limited by 
electron-phonon and electron-electron interactions, which will modify the resonant energies 
and the resonance widths. An assessment of these factors is an important direction for further 
work. The abrupt activation of the waterwheel effect is a strong candidate mechanism behind 
anomalous heating in point contacts [23,24], and can furnish a bias-controlled "switch" for 
current-driven atomic- scale motors. 

The non-equilibrium dynamical response matrix is a useful simple probe into non- 
conservative effects, and is a generalization of the usual description of harmonic motion at 
equilibrium. A complete analysis of the eigenmodes under current requires the inclusion 
of velocity-dependent forces, which can be done perturbatively in the steady state [3,4,9- 
11]. The simpler picture obtained by neglecting the velocity-dependent forces then amounts 
to taking the adiabatic limit of large atomic mass, or small atomic velocities. The eigenmode 
analysis now does not contain, for example, the electronic friction present in our non-adiabatic 
molecular dynamics simulations, but captures rigorously and exactly the physics of steady- 
state conduction in the Born-Oppenheimer limit. 
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